A PRIMER ON EULERIAN COMPUTATIONAL FLUID 
DYNAMICS FOR ASTROPHYSICS 
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ABSTRACT 

We present a pedagogical review of some of the methods employed in Eulerian computational fluid dy- 
namics (CFD). Fluid mechanics is governed by the Euler equations, which are conservation laws for mass, 
momentum, and energy. The standard approach to Eulerian CFD is to divide space into finite volumes 
or cells and store the cell-averaged values of conserved hydro quantities. The integral Euler equations 
are then solved by computing the flux of the mass, momentum, and energy across cell boundaries. We 
review both first-order and second-order flux assignment schemes. All linear schemes are either dispersive 
or diffusive. The nonlinear, second-order accurate total variation diminishing (TVD) approach provides 
high resolution capturing of shocks and prevents unphysical oscillations. We review the relaxing TVD 
scheme, a simple and robust method to solve systems of conservation laws like the Euler equations. 

A 3-D relaxing TVD code is applied to the Sedov- Taylor blast wave test. The propagation of the 
blast wave is accurately captured and the shock front is sharply resolved. We apply a 3-D self-gravitating 
hydro code to simulating the formation of blue straggler stars through stellar mergers and present some 
numerical results. A sample 3-D relaxing TVD code is provided in the appendix. 



Subject headings: hydrodynamics-methods: numerical 

1. Introduction 

Astrophysical structure formation and the dynam- 
ics of astrophysical systems involve nonlinear gas dy- 
namical processes which cannot be modeled analyti- 
cally but require numerical methods. One would like 
to address the challenging problem of star formation 
and how this process produces planetary systems. Ob- 
servations of the X-ray emission from hot gas in galaxy 
clusters, the Sunyaev-Zeldovich effect in the CMB 
spectrum, and the Lyman alpha forest in the spec- 
tra of quasars are only meaningful if we understand 
the gas dynamical processes involved. The evolution 
of complex systems is best modeled using numerical 
simulations. 



A large class of astrophysical problems involve col- 
lisional systems where the mean free path is much 
smaller than all length scales of interest. Hence, one 
can appropriately adopt an ideal fluid description of 
matter where the thermodynamical properties of the 
fluid obey well known equations of state. Conservation 
of mass, momentum, and energy allows one to write 
down the Euler equations which govern fluid mechan- 
ics (See Landau & Lifshitz 1987). This formalism is 
an ideal basis for simulating astrophysical fluids. 

Hydrodynamical simulations are faced with chal- 
lenging problems, but advancements in the field have 
made it an important tool for theoretical astrophysics. 
One of the main challenges in simulating complex fluid 
flows is the capturing of strong shocks, which fre- 
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qucntly occur and play an important role in gas dy- 
namics. The differential Euler equations are ill-defined 
at shock discontinuities where derivatives are infinite. 
Much effort has been devoted to solving this problem 
and a field of work has resulted from it. Computa- 
tional fluid dynamics (CFD) is a powerful approach to 
simulating fluid flow with emphasis on high resolution 
capturing of shocks and prevention of numerical insta- 
bilities. Both Eulerian and Lagrangian methods have 
been developed. 

Lagrangian methods based on smoothed particle 
hydrodynamics (SPH; Gingold & Monaghan 1977; 
Lucy 1977) consider a Monte-Carlo approximation to 
solving the fluid equations, somewhat analogous to TV- 
body methods for the Vlasov equation. SPH schemes 
follow the trajectories of particles of fixed mass which 
represent fluid elements. The Lagrangian forms of the 
Euler equations are solved to determine smoothed fluid 
variables like density, velocity, and temperature. The 
particle formulation does not naturally capture shocks 
and artificial viscosity is added to prevent unphysical 
oscillations. However, the addition of artificial viscos- 
ity broadens shocks over several smoothing lengths and 
degrades the resolution. The Lagrangian approach has 
a large dynamic range in length but not in mass. It 
achieves good spatial resolution in high density regions 
but does poorly in low density regions. SPH schemes 
must smooth over a large number of neighbouring par- 
ticles, making it computationally expensive and chal- 
lenging to implement in parallel. 

The standard approach to Eulerian methods is to 
discrctizc the problem and solve the integral Euler 
equations on a Cartesian grid by computing the flux of 
mass, momentum, and energy across grid cell bound- 
aries. In conservative schemes, the flux taken out of 
one cell is added to the neighbouring cell and this en- 
sures the correct shock propagation. Flux assignment 
schemes based on the total variation diminishing con- 
dition (Hartcn 1983) have been designed for high or- 
der accuracy and high resolution capturing of shocks, 
while preventing unphysical oscillations. The Eule- 
rian approach has a large dynamic range in mass but 
not in length, opposite to that of Lagrangian schemes. 
In general, Eulerian algorithms are computationally 
faster by several orders of magnitude. They are also 
easy to implement and to parallelize. 

The purpose of this paper is to present a pedagog- 
ical review of some of the methods employed in Eu- 
lerian computational fluid dynamics. In §2 we briefly 
review the Euler equations and discuss the standard 



approach to discretizing conservation laws. We de- 
scribe traditional central differencing methods such as 
the Lax-Wendroff scheme in §3 and more modern flux 
assignment methods like the TVD scheme in §4. In 
§5 we review the relaxing TVD method for systems of 
conservation laws like the Euler equations, which has 
been successfully implemented for simulating cosmo- 
logical astrophysical fluids by Pen (1998). In §6 we 
apply a self-gravitating hydro code to simulating the 
formation of blue straggler stars through stellar merg- 
ers. A sample 3-D relaxing TVD code is provided in 
the appendix. 

2. Eulerian Hydrodynamics 

The Euler equations which govern hydrodynamics 
are a system of conservation laws for mass, momen- 
tum, and energy. In differential conservation form, the 
continuity equation, momentum equation, and energy 
equation are given as: 



dp 



d 



de d ,, „. . 
+ — [{ e + P) Vj ]=0 



dt 



(1) 



(2) 



(3) 



We have omitted gravitational and other source terms 
like heating and cooling. The physical state of the fluid 
is specified by its density p, velocity field v, and total 
energy density, 



e + £ 



(4) 



In practice, the thermal energy e is evaluated by sub- 
tracting the kinetic energy from the total energy. For 
an ideal gas, the pressure P(e) is related to the thermal 
energy by the equation of state, 



P=(7-l)e, 



(5) 



where 7 is the ratio of specific heats. Another thermo- 
dynamic variable which is of importance is the sound 
speed c s which is given by 



2 = dp _jp_ 
Cs ~ dp P 



(6) 



The thermodynamical properties of an ideal gas obey 
well known equations of state, which we do not fully 
list here. 
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The differential Euler equations require differen- 
tiable solutions and therefore, are ill-defined at jump 
discontinuities where derivatives are infinite. In the 
literature, nondiffcrentiablc solutions are called weak 
solutions. The differential form gives a complete de- 
scription of the flow in smooth regions, but the integral 
form is needed to properly describe shock discontinu- 
ities. In integral conservation form, the rate of change 
in mass, momentum, and energy is equal to the net 
flux of those conserved quantities through the surface 
enclosing a control volume. For simplicity of notation, 
we will continue to express conservation laws in differ- 
ential form, as a shorthand for the integral form. 

2.1. Computational Fluid Dynamics 

The standard approach to Eulerian computational 
fluid dynamics is to discretize time into discrete steps 
and space into finite volumes or cells, where the con- 
served quantities are stored. In the simplest case, the 
integral Euler equations are solved on a Cartesian cu- 
bical lattice by computing the flux of mass, momen- 
tum, and energy across cell boundaries in discrete time 
steps. Consider the Euler equations in vector differen- 
tial conservation form, 



du 
~dt 



dFjju) 



= 



(7) 



where u = (p, pv Xl pv y , pv z , e) contains the conserved 
physical quantities and F(u) represents the flux terms. 
In practice, the conserved cell-averaged quantities 
u n = u(x n ) and fluxes F n are defined at integer 
grid cell centres x n . The challenge is to use the cell- 
averaged values to determine the fluxes F n+ i/ 2 at cell 
boundaries. 

In the following sections, we describe flux assign- 
ment methods designed to solve conservation laws like 
the Euler equations. For ease of illustration, we begin 
by considering a I-D scalar conservation law, 



du | dF(u) = Q ^ 
dt dx 



(8) 



where F(u) — vu and v is a constant advection veloc- 
ity. Equation (8) is referred to as a linear advection 
equation and has the analytical solution, 



u(x, t) — u(x — vt, 0) . 



(9) 



The linear advection equation describes the transport 
of the quantity u at a constant velocity v. 



In integral flux conservation form, the 1-D scalar 
conservation law can be written as 



d_ 

dt 



rx 2 

•J X\ 



t)dx 



dF{u) 
dx 



dx = 



(10) 



where X\ = x n _ 1 / 2 and x 2 = x n+1 / 2 for our control 
cells. Let F^ +1 , 2 denote the flux of u through cell 
boundary x n+1 / 2 



at time t. Note then that the second 
Ft 1/2 . The rate 



integral is simply equal to , 2 
of change in the cell-integrated quantity / udx is equal 
to the net flux of u through the control cell. For a 
discrete time step, the discretized solution for the cell- 
averaged quantity u n is given by 



ut +At = ui 



r n+l/2 



r n-l/2 



Ax 



At . 



(11) 



The physical quantity u is conserved since the flux 
taken out of one cell is added to the neighbouring cell 
which shares the same boundary. Note that Equa- 
tion (11) has the appearance of being a finite differ- 
ence scheme for solving the differential form of the I- 
D scalar conservation law. This is why the differential 
form can be used as a shorthand for the integral form. 

3. Centered Finite-Difference Methods 

Central-space finite-difference methods have ease of 
implementation but at the cost of lower accuracy and 
stability. For illustrative purposes, we start with a 
simple first-order centered scheme to solve the linear 
advection equation. The discretized solution is given 
by Equation (11) where the fluxes at cell boundaries, 



^n+l/2 



(12) 



are obtained by taking an average of cell-centered 
fluxes F^ = vu l n . The discretized first-order centered 
scheme can be equivalently written as 



ut +At =ul 



F n-1 



2Ax 



At 



(13) 



In this form, the discretization has the appearance 
of using a central difference scheme to approximate 
spatial derivatives. Hence, centered schemes are of- 
ten referred to as central difference schemes. In prac- 
tice when using centered schemes, the discretization is 
done on the differential conservation equation rather 
than the integral equation. 
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This simple scheme is numerically unstable and we 
can show this using the von Neumann linear stability 
analysis. Consider writing u{x,t) as a discrete Fourier 
series: 



u„ = 



N/2 

N ^ 

k=-N/2 



f ( 2irikn\ 

Ckexp \~nr) 



(14) 



where N is the number of cells in our periodic box. In 
plane-wave solution form, we can write this as 



N/2 

N ^ 



c k cxp 



k=-N/2 



2m(kn — tot) 
N 



(15) 



where c£ are the Fourier series coefficients for the ini- 
tial conditions u(x,0). Equivalcntly, the time evolu- 
tion of the Fourier series coefficients in Equation (14) 
can be cast into a plane-wave solution of the form, 
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cxp 



-2iriujt\ 



N / 



(16) 



where the numerical dispersion relation u)(k) is com- 
plex in general. The imaginary part of u> represents the 
growth or decay of the Fourier modes while the real 
part describes the oscillations. A numerical scheme 
is linearly stable if Im(w) < 0. Otherwise, the Fourier 
modes will grow exponentially in time and the solution 
will blow up. 

The exact solution to the linear advection equation 
can be expressed in the form of Equation (9) or by 
a plane-wave solution where the dispersion relation is 
given by uj — vk. The waves all travel at the same 
phase velocity uo /k = v in the exact case. 

The centrally discretized linear advection equation 
(Equation 13) is exactly solvable. After m times steps, 
the time evolution of the independent Fourier modes 
is given by 



(17) 



where A = vAt/Ax and 
sion relation is given by 



2nkAx/N. The disper- 



N 
2irAt 



tan ^Asin^) + -ln(l + A 2 sin 2 



(18) 

For any time step At > 0, the imaginary part of to will 
be > 0. The Fourier modes will grow exponentially in 
time and the solution will blow up. Hence, the first- 
order centered scheme is numerically unstable. 



3.1. Lax-Wendroff Scheme 

The Lax-Wendroff scheme (Lax & Wendroff 1960) is 
second-order accurate in time and space and the idea 
behind it is to stabilize the unstable first-order scheme 
from the previous section. Consider a Taylor series 
expansion for u(x, t + At): 

u(x, t + At) — u(x, t) 

du . d 2 u At 2 m . A 
+ — At + — — + 0{At 3 ) , (19) 



dt 



dt 2 2 



and replace the time derivatives with spatial deriva- 
tives using the conservation law to obtain 

u(x, t + At) — u(x, t) 

dF d fdFdFdu\At 2 ^ /A ,, 
-^ M+ lTx{^d-x)— +0iAt) - (20) 

For the linear advection equation, the eigenvalue of 
the flux Jacobian is dF/du = v. Discretization using 
central differences gives 



,t+At „ „t r n+l r n-l 



= u„- 



At 



+ 



2Ax 

Ax Ax I 2Ax 



■ (21) 



In conservation form, the solution is given by Equation 
(11), where the fluxes at cell boundaries are defined as 



1 



j At 



F n+l/2 - 2 + F n) ~ ( F n+1 ~ F n) ■ ( 22 ) 

Compare this with the boundary fluxes for the first- 
order scheme (Equation 12). The Lax-Wendroff 
scheme obtains second-order fluxes, 



r(2) _ 



dF dF At 

du dx 2 



(23) 



by modifying the first-order fluxes with a second- 
order correction. 

The stability of the Lax-Wendroff scheme to solve 
the linear advection equation can also be determined 
using the von Neumann analysis. The discretized Lax- 
Wendroff equation (Equation 21) is exactly solvable 
and after m time steps, the Fourier modes evolve ac- 
cording to 

cf'= [1 - A 2 (l -cos0) -i\sm(j)] m cl , (24) 
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Fig. 1. — The phase error Re(Aw) and the amplica- 
tion factor Im(Aw) for the Lax-Wcndroff scheme with 
parameters N — 100, v = 1, and A = 0.9. 

where A = vAt/ Ax is called the Gourant number and 
<t> = 2-nkAx/N. The dispersion relation is given by 



N 
2-kAI 



tan 



Asin^ 



iN , 

H t— In 

4nAt 



1 - A 2 (l -cos</>) 



1-4A 2 (1-A 2 )sin 4 



(25) 



It is important to note three things. First, the Lax- 
Wendroff scheme is conditionally stable provided that 
Im(w) < 0, which is satisfied if 



vAt 
Ax 



< 1 



(26) 



This constraint is a particular example of a general 
stability constraint known as the Courant-Friedrichs- 
Lewyor (CFL) condition. The Courant number A is 
also referred to as the CFL number. Second, for A = 1 
the dispersion relation is exactly identical to that of 
the exact solution and the numerical advection is ex- 
act. This is a special case, however, and it does not test 
the ability of the Lax-Wcndroff scheme to solve general 
scalar conservation laws. Normally, one chooses A < 1 
to satisfy the CFL condition. Lastly, for A < 1 the dis- 
persion relation uu{k) for the Lax-Wendroff solution is 



Fig. 2. — Lax-Wendroff scheme used to linearly advect 
a square wave (solid line) once (dashed line) and ten 
times (dotted line) through a box of 100 grid cells at 
speed v = 1. 

different from the exact solution where u a — vk. The 
dispersion relation relative to the exact solution can 
be parametrized by 



AiO = LO — L0 o 



(27) 



The second-order truncation of the Taylor series 
(Equation 20) results in a phase error Re(Aw) which is 
a function of frequency. In the Lax-Wcndroff solution, 
the waves are damped and travel at different speeds. 
Hence the scheme is both diffusive and dispersive. 

In Figure (1) we plot the phase error Re(Aw) and 
the amplification term Im(Aw) for the Lax-Wcndroff 
scheme with parameters N = 100, v = 1, and A = 0.9. 
A negative value of Re(Aw) represents a lagging phase 
error while a positive value indicates a leading phase 
error. For the chosen CFL number, the high frequency 
modes have the largest phase errors but they are highly 
damped. Some of the modes having lagging phase er- 
rors are not highly damped. We will subsequently see 
how this becomes important. 

A rigourous test of the 1-D Lax-Wendroff scheme 
and other flux assignment schemes we will discuss is 
the linear advection of a square wave. The challenge is 
to accurately advect this discontinuous function where 
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Fig. 3. — The phase error Re(Aw) and the amplication 
term Im(Aw) for the Lax-Wendroff scheme (boxes) 
and the first-order upwind scheme (crosses) with pa- 
rameters N = 100, v = 1, and A = 0.9. 



Fig. 4. — First-order upwind scheme used to linearly 
advect a square wave (solid line) once (dashed line) 
and ten times (dotted line) through a box of 100 grid 
cells at speed v = 1. 



the edges mimic Riemann shock fronts. In Figure (2) 
we show how the Lax-Wendroff scheme does at advect- 
ing the square wave once (dashed line) and ten times 
(dotted line) through a periodic box of 100 grid cells 
at speed v = 1 and A = 0.9. Note that this scheme 
produces numerical oscillations. Recall that a square 
wave can be represented by a sum of Fourier or sine 
waves. These waves will be damped and disperse when 
advected using the Lax-Wendroff scheme. Figure (1) 
shows that the modes having lagging phase errors are 
not damped away. Hence, the Lax-Wendroff scheme is 
highly dispersive and the oscillations in Figure (2) are 
due to dispersion. We leave it as an exercise for the 
reader to advect a sine wave using the Lax-Wendroff 
scheme. Since there is only one frequency mode in this 
case, there will be no spurious oscillations due to dis- 
persion, but a phase error will be present. For a com- 
prehensive discussion on the family of Lax-Wendroff 
schemes and other centered schemes, see Hirsch (1990) 
and Laney (1998). 



4. Upwind Methods 

Upwind methods take into account the physical na- 
ture of the flow when assigning fluxes for the discrete 
solution. This class of flux assignment schemes, whose 
origin dates back to the work of Courant, Isaason, & 
Reeves (1952), has been shown to be excellent at cap- 
turing shocks and also being highly stable. 

We start with a simple first-order upwind scheme 
to solve the linear advection equation. Consider the 
case where the advection velocity is positive and flow 
is to the right. The flux of the physical quantity u 
through the cell boundary x n+ i/ 2 will originate from 
cell n. The upwind scheme proposes that, to first- 
order, the fluxes F^ +1 ^ 2 at cell boundaries be taken 
from the cell-centered fluxes F* = vu l n , which is in the 
upwind direction. If the advection velocity is negative 
and flow is to the left, the boundary fluxes i^* +1 / 2 
are taken from the cell-centered fluxes = vu^+i- 

The first-order upwind flux assignment scheme can be 
summarized as follows: 



r n+l/2 



F* if v > , 



(28) 
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Unlike central difference schemes, upwind schemes are 
explicitly asymmetric. 

The CFL condition for the first-order upwind 
scheme can be determined from the von Neumann 
analysis. We consider the case of a positive advection 
velocity. After m time steps, the Fourier modes evolve 
according to 



„mAt 



[1 — A(l — cos0) — iAsin< 



(29) 



where A = vAt/Ax and = 2irkAx/N. The disper- 
sion relation is given by 



U) = 



N 
2nAt 



tan 



A sin ( 



1 — A(l — cos</>) 



+ 



iN 
AirAt 



hi 



1 - 4A(1 - A) sin 2 £ 



(30) 



The CFL condition for solving the linear advection 
equation with this scheme is to have A < 1, identical to 
that for the Lax-Wendroff scheme. For A < 1 the dis- 
persion relation u(k) for the first-order upwind scheme 
is different from the exact solution where u = vk. 
This scheme is both diffusive and dispersive. Since it 
is only first-order accurate, the amount of diffusion is 
large. In Figure (3) we compare the dispersion rela- 
tion of the upwind scheme to that of the Lax-Wendroff 
scheme. The Fourier modes in the upwind scheme also 
have phase errors but they will be damped away. The 
low frequency modes which contribute to the oscilla- 
tions in the Lax-Wendroff solution are more damped 
in the upwind solution. Hence, one does not expect to 
see oscillations resulting from phase errors. 

In Figure (4) we show how the first-order upwind 
scheme does at advecting the Riemann shock wave. 
This scheme is well-behaved and produces no spurious 
oscillations, but since it is only first-order, it is highly 
diffusive. The first-order upwind scheme has the prop- 
erty of having monotonicity preservation. When ap- 
plied to the linear advection equation, it does not allow 
the creation of new extrema in the form of spurious os- 
cillations. The Lax-Wendroff scheme does not have the 
property of having monotonicity preservation. 

The flux assignment schemes that we have discussed 
so far are all linear schemes. Godunov (1959) showed 
that all linear schemes are either diffusive or dispersive 
or a combination of both. The Lax-Wendroff scheme is 
highly dispersive while the first-order upwind scheme 
is highly diffusive. Godunov's theorem also states 
that linear monotonicity preserving schemes are only 
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Fig. 5. — The TVD scheme using the minmod flux 
limiter is applied to the advection of a square wave. 

first-order accurate. In order to obtain higher order 
accuracy and prevent spurious oscillations, nonlinear 
schemes are needed to solve conservation laws. 

4.1. Total Variation Diminishing Schemes 

Harten (1983) proposed the total variation di- 
minishing (TVD) condition which guarantees that 
a scheme have monotonicity preservation. Applying 
Godunov's theorem, we know that all linear TVD 
schemes are only first-order accurate. In fact, the only 
linear TVD schemes are the class of first-order up- 
wind schemes. Therefore, higher order accurate TVD 
schemes must be nonlinear. 

The TVD condition is a nonlinear stability condi- 
tion. The total variation of a discrete solution, defined 
as 



TV(u*) 



N 

1=1 



+1 



(31) 



is a measure of the overall amount oscillations in u. 
The direct connection between the total variation and 
the overall amount of oscillations can be seen in the 
equivalent definition 

TV(u') = 2 
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Fig. 6. — The TVD scheme using the superbee flux 
limiter is applied to the advection of a square wave. 



Fig. 7. — The TVD scheme using the Van Leer flux 
limiter is applied to the advection of a square wave. 



where each maxima is counted positively twice and 
each minima counted negatively twice (See Laney 
1998). The formation of spurious oscillations will con- 
tribute new maxima and minima and the total varia- 
tion will increase. A flux assignment scheme is said to 
be TVD if 

TV(u t+At ) < TV{u l ) , (33) 

which signifies that the overall amount of oscillations is 
bounded. In linear flux-assignment schemes, the von 
Neumann linear stability condition requires that the 
Fourier modes remain bounded. In nonlinear schemes, 
the TVD stability condition requires that the total 
variation diminishes. 

We now describe a nonlinear second-order accurate 
TVD scheme which builds upon the first-order mono- 
tone upwind scheme described in the previous section. 
The second-order accurate fluxes ^* +1 / 2 a * ceu bound- 
aries are obtained by taking first-order fluxes 
from the upwind scheme and modifying it with a sec- 
ond order correction. First consider the case where the 
advection velocity is positive. The first-order upwind 
flux F^y 2 comes from the averaged flux F„ in cell n. 



We can define two second-order flux corrections, 

F t 



F n-1 



F n+1 



(34) 



(35) 



using three local cell-centered fluxes. We use cell n 
and the cells immediately left and right of it. If the 
advection velocity is negative, the first-order upwind 
flux comes from the averaged flux F^ +1 in cell n + 1. 
In this case, the second-order flux corrections, 



F n+1 



F n+2 



F n+1 



(36) 



(37) 



are based on cell n+1 and the cells directly adjacent to 
it. Near extrema where the corrections have opposite 
signs, we impose no second-order correction and the 
flux assignment scheme reduces to first-order. A flux 
limiter <f> is then used to determine the appropriate 
second-order correction, 



^+\ /2 ,AF^ 1/2 ), 



(38) 
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which still maintains the TVD condition. The second- 
order correction is added to the first-order fluxes to get 
second-order fluxes. The first-order upwind scheme 
and second-order TVD scheme will be referred to 
as monotone upwind schemes for conservation laws 
(MUSCL). 

Time integration is performed using a second-order 
Rungc-Kutta scheme. We first do a half time step, 




(39) 



using the first-order upwind scheme to obtain the half- 
step values u t+At / 2 . A full time step is then computed, 



ul +At = ut- 



n+1/2 



-^t+At/2 
1/2 



Ax 



At , (40) 



using the TVD scheme on the half-step fluxes ^*+ At / 2 
The reader is encouraged to show that is second-order 
accurate. 

We briefly discuss three TVD limiters. The min- 
mod flux limiter chooses the smallest absolute value 
between the left and right corrections: 

minmod(a,b) = ^[sign(a) + sign(6)] min(|a|, |6|) . 

(41) 

The superbee limiter (Roe 1985) chooses between the 
larger correction and two times the smaller correction, 
whichever is smaller in magnitude: 

minmod(a, 2b) if \a\ > \b\ , 
superbee(a, b) = ^ (42) 
minmod(2a, b) otherwise . 

The Van Leer limiter (Van Leer 1974) takes the har- 
monic mean of the left and right corrections: 



vanleer(a, b) 



2ab 
a + b ' 



(43) 



The minmod limiter is the most moderate of all 
second-order TVD limiters. In Figure (5) we see that 
it does not do much better than first-order upwind 
for the square wave advection test. Superbee chooses 
the maximum correction allowed under the TVD con- 
straint. R is especially suited for piece- wise linear 
conditions and is the least diffusive for this particular 
test, as can be seen in Figure (6). Note that no ad- 
ditional diffusion can be seen by advecting the square 
wave more than once through the box. It can be shown 



that the minmod and superbee limiters are extreme 
cases which bound all other second-order TVD lim- 
iters. The Van Leer limiter differs from the previous 
two in that it is analytic. This symmetrical approach 
falls somewhere inbetween the other two limiters in 
terms of moderation and diffusion, as can be seen in 
Figure (7). It can be shown that the CFL condition 
for the second-order TVD scheme is to have A < 1. 
For a comprehensive discussion on TVD limiters, see 
Hirsch (1990) and Laney (1998). 

5. Relaxing TVD 

We now describe a simple and robust method to 
solve the Euler equations using the monotone upwind 
scheme for conservation laws (MUSCL) from the pre- 
vious section. The relaxing TVD method (Jin & Xin 
1995) provides high resolution capturing of shocks us- 
ing computationally inexpensive algorithms which are 
straightforward to implement and to parallelize. It 
has been successfully implemented for simulating cos- 
mological astrophysical fluids by Pen (1998). 

The MUSCL scheme is straightforward to apply to 
conservation laws like the advection equation since the 
velocity alone can be used as a marker of the direc- 
tion of flow. However, applying the MUSCL scheme to 
solve the Euler equations is made difficult by the fact 
that the momentum and energy fluxes depend on the 
pressure. In order to determine the direction upwind of 
the flow, it becomes necessary to calculate the flux Ja- 
cobian eigenvectors using Ricmann solvers. This step 
requires computationally expensive algorithms. The 
relaxing TVD method offers an attractive alternative. 

5.1. 1-D Scalar Conservation Law 

We first present a motivation for the relaxing 
scheme by again considering the 1-D scalar conser- 
vation law. The MUSCL scheme for solving the linear 
advection equation is explicitly asymmetric in that it 
depends on the sign of the advection velocity. We now 
describe a symmetrical approach which applies to a 
general advection velocity. 

The flow can be considered as a sum of a right- 
moving wave u R and a left-moving wave u L . For a 
positive advection velocity, the amplitude of the left- 
moving wave is zero and for a negative advection ve- 
locity, the amplitude of the right-moving wave is zero. 
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In compact notation, the waves can be denned as: 



u L = 



l + v/c 



1 — v/c 



u , 



(44) 



(45) 



where c — \v\. The two waves are traveling in opposite 
directions with advection speed c and can be described 
by the advection equations: 



du R d , p . 



dt 



, 



(46) 



d 



(cu L ) = 



(47) 



du L 

dt dx 1 

The MUSCL scheme is straightforward to apply to 
solve Equations (46) and (47) since the upwind direc- 
tion is left for the right-moving wave and right for the 
left-moving wave. The 1-D relaxing advection equa- 
tion then becomes 



du dF R dF L 
dt dx dx 



(48) 



where F R = 



,R 



and F L = 



For the discretized 



solution given by Equation (11), the boundary fluxes 

^n+l/2 are n0W a SUnl °^ tne nuxes F„+l/2 an< ^ 

from the right-moving and left-moving waves, respec- 
tively. Note that the relaxing advection equation will 
correctly reduce to the linear advection equation for 
any general advection velocity. 

Using this symmetrical approach, a general algo- 
rithm can be written to solve the linear advection 
equation for an arbitrary advection velocity. This 
scheme is indeed inefficient for solving the linear ad- 
vection equation since one wave will have zero am- 
plitude. However, the Euler equations can have both 
right-moving and left-moving waves with non-zero am- 
plitudes. 

5.2. 1-D Systems of Conservation Laws 

We now discuss the 1-D relaxing TVD scheme and 
later generalize it to higher spatial dimensions. Con- 
sider a 1-D system of conservation laws, 



du 
~dt 



dF{u) 

dx 



, 



(49) 



where for the Euler equations, we have u = (p,pv,e) 
and F(u) the corresponding flux terms. We now re- 
place the vector conservation law with the relaxation 



system 



du d , . 

m + d-x {cw) 

dw d . . 
-dt + d~x {cU) 



, 



(50) 



(51) 



where c(x, t) is a free positive function called the freez- 
ing speed. The relaxation system contains two coupled 
vector linear advection equations. In practice, we set 
w = F(u)/c and use it as an auxiliary vector to cal- 
culate fluxes. Equation (50) reduces to our 1-D vector 
conservation law and Equation (51) is a vector conser- 
vation law for w. 

In order to solve the relaxed system, we decouple 
the equations through a change of variables: 



u 



u + w 



u — w 



u = 



which then gives us 

du R 
~df 

du L 
~df 



d~x {cU) 







(52) 



(53) 



(54) 



(55) 



Equations (54) and (55) are vector linear advection 
equations, which can be interpreted as right-moving 
and left-moving flows with advection speed c. Note 
the similarity with their scalar counterparts, Equa- 
tions (46) and (47). The 1-D vector relaxing conser- 
vation law for u becomes 



dF l 



du dF 1 

dt dx dx 







(56) 



where F R = cu R and F L = cu L . The vector relaxing 
equation can now be solved by applying the MUSCL 
scheme to Equations (54) and (55). Again, note the 
similarity between the vector relaxing equation and its 
scalar counterpart, Equation (48). 

The relaxed scheme is TVD under the constraint 
that the freezing speed c be greater than the charac- 
teristic speed given by the largest eigenvalue of the flux 
Jacobian dF(u)/du. For the Euler equations, this is 
satisfied for 

c= \v\ +c s . (57) 

Jin & Xin (1995) considered the freezing speed to be 
a positive constant in their relaxing scheme while we 
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generalize it to be a positive function. Time integra- 
tion is again performed using a second-order Runge- 
Kutta scheme and the time step is determined by sat- 
isfying the CFL condition, 



Ax 



< 1 . 



(58) 



Note that a new freezing speed is computed for each 
partial step in the Runge-Kutta scheme. The CFL 
number A = c max At/ Ax should be chosen such that 
Cmax will be larger than max(c^) and max(cn +At//2 ). 

We now summarize the steps needed to numerically 
solve the 1-d Euler equations. At the beginning of 
each partial step in the Runge-Kutta time integra- 
tion scheme, we need to calculate the cell-averaged 
variables defined at grid cell centres. First for the 
half time step, we calculate the fluxes -F( u ri) and the 
freezing speed c n . We then set the auxiliary vector 
w n = F(v,n)/Cn and construct the right-moving waves 
u R ' 1 and left-moving waves u^'*. The half time step 
is given by 



u 



t+At/2 



Ul, 



n+1/2 



— F 



-1/2 



Ax 



where 

F n+l/2 = F n+1/2 ~ F n+1/2 ' ( 60 ) 

The first-order upwind scheme is used to compute 
fluxes at cell boundaries for the right-moving and left- 
moving waves. For the full time step, we construct the 
right-moving waves u R ' t+At ^ 2 and left-moving waves 
u L,t+At/2^ ugm g j. ne half-step values of the appropri- 
ate variables. The full time step, 



ut +At = ui 



? t+At/2 
n+1/2 



p t+At/2' 
r n-1/2 



Ax 



At 



(61) 



is computed using the second-order TVD scheme. This 
completes the updating of it* to u t+At . 

We have found that a minor modification to the 
implementation described above gives more accurate 
results. Consider writing the flux of the right-moving 
and left-moving waves as: 



cG 1 



F = cG L 



(62) 
(63) 



where G R is the flux of fi R = u R /c and G L is the 
flux of fi L = u L /c. The linear advection equations for 



fi R and fi L are similar to Equations (54) and (55), but 
where we replace u R with fi R and u L with fi L . For 
each partial step in the Runge-Kutta scheme, the net 
fluxes at cell boundaries are then taken to be 



? n+l/2 — c n+l/2(G R +l/2 



Gn+1/2) 



(64) 



where we use c n +i/2 = {c n +i + c„)/2. In practice, this 
modified implementation has been found to resolve 
shocks with better accuracy in certain cases. Note 
that the two different implementations of the relaxing 
TVD scheme are identical when a constant freezing 
speed is used. 

5.3. Multi-Dimensional Systems of Conserva- 
tion Laws 

The 1-D relaxing TVD scheme can be generalized 
to higher dimensions using the dimensional splitting 
technique by Strang (1968). In three dimensions, the 
Euler equations can be dimensionally split into three 
separate 1-D equations which are solved sequentially. 
Let the operator Li represent the updating of u l to 
U t+At ky j nc i uc jing the flux in the i direction. We first 
complete a forward sweep, 



,t+At 



— L Z LyL 



and then perform a reverse sweep 



t+2At 



L x L y L z u 



t+At 



(65) 



(66) 



using the same time step At to obtain second-order 
accuracy. We will refer to the combination of the for- 
ward and reverse sweeps as a double sweep. 

A more symmetrical sweeping pattern can be used 
by permutating the sweeping order when completing 
the next double time step. The dimensional splitting 
or operator splitting technique can be summarized as 
follows: 



u t 2=u t 1+ 2A tl _ 



— L x L y L z L z LyL x u 



ti 



u t 3 = u t 2 +2At 2 = LzLxLyLyLxLzU t2 



i 3 ~^~ 3 — LyL z L x L x L z LyU 3 , 



(67) 
(68) 
(69) 



where At\, At2, and At3 are newly determined time 
steps after completing each double sweep. 

The CFL condition for the 3-D relaxing TVD 
scheme is similarly given by Equation (58), but with 



Cmax — ^^[(^xjmaxj (Cy)max? (^z)r 



(70) 
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Fig. 8. — Sedov- Taylor blast wave test conducted in a 
box with 256 3 cells. The data points are taken from 
a random subset of cells and the solid lines are the 
analytical self-similar solutions. 

where c, = \vi\ + c s . Note that since max(|wi|) is on 
average a factor of \/3 smaller than max(|v|), a di- 
mensionally split scheme can use a longer time step 
compared to an un-split scheme. 

The dimensional splitting technique also has other 
advantages. The decomposition into a 1-D problem al- 
lows one to write short 1-D algorithms, which are easy 
to optimize to be cache efficient. A 3-D hydro code is 
straightforward to implement in parallel. When sweep- 
ing in the x direction, for example, one can break up 
the data into 1-D columns and operate on the indepen- 
dent columns in parallel. A sample 3-D relaxing TVD 
code, implemented in parallel using OpenMP direc- 
tives, is provided in the appendix. 

6. Sedov- Taylor Blast Wave Test for 3-D Hy- 
dro 

A rigourous and challenging test for any 3-D Eule- 
rian or Lagrangian hydrodynamic code is the Sedov- 
Taylor blast wave test. We set up the simulation box 
with a homogeneous medium of density pi and neg- 
ligible pressure and introduce a point-like supply of 
thermal energy E Q in the centre of the box at time 



Fig. 9. — A closeup of the Sedov- Taylor blast wave. 
The resolution of the shock front is roughly two grid 
cells and the anisotropic scatter is less than one grid 
cell. 

t = 0. The challenge is to accurately capture the 
strong spherical shock wave which sweeps along ma- 
terial as it propagates out into the ambient medium. 
The Sedov- Taylor test is used to model nuclear-type 
explosions. In astrophysics, it is often used as a basic 
setup to model supernova explosions and the evolution 
of supernova remnants (See Shu 1992). 

The analytical Sedov solution uses the self-similar 
nature of the blast wave expansion (See Landau & Lif- 
shitz 1987). Consider a frame fixed relative to the cen- 
tre of the explosion. The spherical shock front propa- 
gates outward and the distance from the origin is given 

by 

/ p 4-2 \ V 5 

r sh (t)=£ (^-) , (71) 

where £ = 1-15 for an ideal gas with 7 = 5/3. The 
velocity of the shock v s h = dr s h / dt is given by 

v sh {t) = \ T -^l . (72) 

Since the ambient medium has negligible pressure, the 
shocks will be very strong. The density P2, velocity v 2} 
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and pressure P2 directly behind the shock front are: 
7+1 



are given as: 



Pi 



V2 



Po 



7-1 



7+1 



7+1 



Pi 



V s h , 



Piv 2 sh 



(73) 



(74) 



(75) 



The immediate post-shock gas density is constant in 
time, while the shocked gas velocity V2 and pressure P2 
decrease as t~ 3 ^ 5 and i~ 6 / 5 , respectively. The full ana- 
lytical Sedov- Taylor solutions can be found in Landau 
& Lifshitz (1987). 

The 3-D relaxing TVD code based on the van Leer 
flux limiter is applied to capturing the Sedov- Taylor 
blast wave. We set up a box with 256 3 cells and con- 
stant initial density p\ — 1. At time t — 0, we inject a 
supply of thermal energy E = 10 5 into one cell at the 
centre of the box. The simulation is stopped at time 
t = 283, in which the shock front has propagated out to 
a distance of r s h = 110 cells from the centre. In Figure 
(8) and (9) we plot the radial distributions of density, 
momentum, and pressure, normalized to P2, P2V2, and 
P2 respectively. The data points are taken from a ran- 
dom subset of cells and the solid lines are the analytical 
Sedov- Taylor solutions. Despite solving a spherically 
symmetric problem on an explicitly non-rotationally 
invariant Cartesian grid, the anisotropic scatter in the 
results is small. The distance of the shock front from 
the centre of the explosion as a function of time is in- 
deed given by Equation (71), demonstrating that the 
3-D relaxing TVD code ensures the correct shock prop- 
agation. The resolution of the shock front is roughly 
two grid cells. The numerical shock jump values of P2, 
V2 , and P2 are resolution dependent and come close to 
the theoretical values for our test with 256 3 cells. We 
leave it as an exercise for the reader to test the code 
using the minmod and superbee flux limiters. 

7. Self- Gravitating Hydro for Astrophysical 
Applications 

For astrophysical applications, both hydrodynami- 
cal and gravitational forces are included. The gravi- 
tational forces arise from the self-gravity of the fluid 
and can also come from an external field. The Euler 
equations with the gravitational source terms included 



(76) 



d(pvi) d 



+ fyTiPViVj + PSij) = -p— , (77) 

where <j> is the gravitational potential. Poisson's equa- 
tion, 

V 2 cj) = 4nGp , (79) 

relates the gravitational potential to the density field. 
The general solution can be written as 

<j>(x) = j p(x')w(x - x')d 3 x' , (80) 

where the kernel is given by 



w(x) = 



G 



x 



(81) 



In the discrete case, the integral in Equation (80) be- 
comes a sum and Poisson's equation can be solved us- 
ing fast Fourier transforms (FFT) to do the convolu- 
tion. The forces are then calculated by finite differenc- 
ing the potential (See Hockney & Eastwood 1988). 

The addition of gravatitational source terms in the 
Euler equations is easily handled using the operator 
splitting technique described in §5.3. Consider the 
double sweep: 



,t+2At 



L x LyL z GGL z LyL x u 



(82) 



where the operator Li represents the updating of u by 
including the flux in the i direction and the operator 
G represents the gravitational acceleration of the fluid. 
During the gravitational step, the flux terms in the 
Euler equations are ignored. The density distribution 
does not change and only the fluid momenta and total 
energy density are updated. 

7.1. Astrophysical Formation of Blue Strag- 
glers Through Stellar Collisions 

The stellar density in the cores of globular and open 
clusters is high enough for stellar collisions to take 
place with significant frequency (Hills & Day 1976). 
Current observations and simulations suggest that the 
merger of two main sequence stars produces a blue 
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straggler (Sills et al. 1997; Sandquist, Bolte, & Hcrn- 
quist 1997). The blue stragglers are out-lying main 
sequence stars which lie beyond the main sequence 
turnoff in the colour-magnitude diagram (CMD) of a 
star cluster. The blue stragglers are more massive, 
brighter, and bluer than the turnoff stars. Since more 
massive stars evolve faster than lower mass stars and 
are not expected to lie beyond the turnoff, this sug- 
gests that blue stragglers must have formed more re- 
cently. 

In principle the merger of two main sequence stars 
can produce a young remnant star provided that sig- 
nificant mixing occurs in the process. The mixing pro- 
duces a higher hydrogen fraction in the core of the 
remnant than that of the parent stars which have al- 
ready burnt most of the hydrogen to helium in their 
cores. Bcnz & Hills (1987) used low resolution SPH 
simulations with ~ 10 3 particles to simulate the merg- 
ing of n = 3/2 polytropes and found that they fully 
mixed. However, medium resolution SPH simulations 
with ~ 10 4 particles of n = 3/2 or n = 3 polytropes 
showed only weak mixing (Lombardi, Rasio, & Shapiro 
1996; Sandquist, Bolte, & Hernquist 1997). It is worth 
noting that n = 3/2 polytropes are more representa- 
tive of low mass main sequence stars with large convec- 
tive envelopes while n = 3 polytropes resemble main 
sequence stars near the turnoff which have little mass 
in their convective envelopes. High resolution SPH 
simulations involving ~ 10 5 — 10 6 particles have now 
been applied to simulating stellar collisions (Sills et a. 
2002). 

The merging stars process is mostly subsonic and 
strong shocks are not expected. In the absence of 
shocks, SPH particles will follow flow lines of constant 
entropy due to the Lagrangian nature of the method. 
As a result, the particles may experience sedimenta- 
tion. In addition, the mixing can also depend on the 
adopted smoothing length and the form of artificial 
viscosity. For a SPH fluid, the Reynolds number is 
of order (Np/Ng) 1 ^ 3 , where N p is the total number of 
particles and N s is the number of particles over which 
the smoothing is done. For N p <~ 10 5 and N s ~ 10 2 , 
the Reynolds number is ~ 10. However, a fluid with a 
low Reynolds number will tend to experience laminar 
flow. Hence, SPH may under mix. 

It is a worthwhile exercise to model the merging 
process using Eulerian hydrodynamical simulations. 
The differences between Eulerian and Lagrangian ap- 
proaches may lead to very different results on mixing. 
As of present, no such work has been reported in the 
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Fig. 10. — The advection of a self-gravitating poly- 
trope in a periodic box with 256 3 cells. We compare 
the mass and entropy profiles of the initial (solid line) 
and advected polytrope after 1000 timesteps in which 
the polytrope has moved 256 cells in each direction. 

literature. 

7.2. Numerical Method 

We consider the off-axis collision of two main se- 
quence stars with M = 0.8 M© and R — 0.955 Rq, 
which arc modeled using n = 3 polytropes. A poly- 
trope with polytropic index n has equilibrium density 
and pressure profiles which are related by 

P oc p 1+1/n . (83) 

The density profile is determined by solving the Lane- 
Emden equation (See Chandrasekhar 1957). We adopt 
an ideal gas equation of state with adiabatic index 
7 = 5/3. Note that for an n = 3 polytrope, 90% 
of the total mass is contained within r < 0.5R. We 
define the dynamical time to be, 



Tdyn — 




where p is the average density. For the chosen parent 
stars, the dynamical time is approximately one physi- 
cal hour. 
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Fig. 11. — Snapshots of the merging process taken at 
time t = 0, 2, 4, and 8rd yn - The density contours are 
spaced logarithmically with 2 per decade and covering 
3 decades down from the maximum. 



The collision is simulated in a box with 512 x 512 x 
256 cells and the orbital plane coincides with the x — y 
plane. Initially each parent star has a radius of 96 grid 
cells. The stars are set up on zero-energy parabolic 
orbits with a periccntre separation equal to 0.25 R. 
The initial trajectories are calculated assuming point 
masses. In an Eulerian simulation, the vacuum cannot 
have zero density. We set the minimum density of the 
cells to be 10~ 8 of the central density of the parent 
stars. The hydrodynamics is done in a non-periodic 
box with vacuum boundary conditions. 

7.3. Numerical Results 

A non-trivial test of a self-gravitating Eulerian hy- 
dro code is the advection of an object in hydrostatic 
equilibrium. The challenge is to maintain the equilib- 
rium profile over a large number of time steps. One of 
the parent stars is placed in a periodic box with 256 3 
cells and given some initial momentum. We make the 
test rigorous by having the polytrope move in all three 
directions. In Figure 10 we compare the mass and en- 
tropy profiles of the initial and advected polytrope. 



The entropic variable A = P/p 1 



sd in place of 



the specific entropy. The parameter A a is defined to 



Fig. 12. — Thermodynamic profiles of the merger rem- 
nant (solid) and the parent stars (dashed). Units are 
cgs. The radial plot is included for comparison. 



be the minimum entropy of the parent polytrope. Af- 
ter 1000 timesteps in which the polytrope has moved 
256 cells in each direction, the advected polytrope has 
still retained its equilibrium profile. Shock heating can 
occur in the outer envelope as the polytrope moves 
through the false vacuum. However, by setting the 
density of the false vacuum to be 10~ 8 of the central 
density of the polytrope, we can minimize the spurious 
shock heating. 

In Figure 1 1 we show four snapshots of the merging 
process taken at time t = 0, 2, 4, and & Tdy n . The 2-D 
density maps are created by averaging over 4 planes 
taken about the orbital mid-plane. The density con- 
tours are spaced logarithmically with 2 per decade and 
covering 3 decades down from the maximum. The par- 
ent stars are initially separated by 3.75 R and placed 
on zero-energy orbits with a pericentre separation of 
0.25i?. During the collision process, the outer en- 
velopes of the parent stars are shock heated and mate- 
rial gets ejected. In less than 10 Td yn , the merger rem- 
nant establishes hydrostatic equilibrium. The merger 
remnant is a rotating oblate with mass approximately 
90% of the combined mass of the parent stars. A large 
fraction of the mass loss is due to the vacuum bound- 
ary conditions. Ejected material do not have the op- 
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portunity to fall back onto the merger remnant. How- 
ever, the additional mass loss in the envelope does not 
present a problem since we are interested in the ques- 
tion of mixing in the interior of the star. 

In Figure 12 we plot the thermodynamic profiles of 
the merger remnant and the parent stars. The central 
density and pressure in the core of the merger remnant 
is lower than the corresponding values in the parent 
stars by approximately half. The entropy floor has 
risen by a factor of 1.6. Shock heating is expected to 
be minimal in the core so a change in entropy suggests 
that some mixing has taken place. However, it is diffi- 
cult to quantify the amount of mixing from examining 
the thermodynamic profiles alone. 

7.4. Future Work 

To help address the question of mixing, we are im- 
plementing a particle-mesh (PM) scheme where test 
particles can be used to track passively advected quan- 
tities such as chemical composition. Initially, each par- 
ent star is assigned a large number of particles with 
known chemical composition. The test particles are 
passively advected along velocity field lines. For each 
time step, the velocity of each particle is interpolated 
from the grid using a "cloud-in-cell" (CIC) scheme 
(Hockncy & Eastwood 1988) and the equations of mo- 
tions are solved using second-order Runge-Kutta in- 
tegration. The CIC interpolation scheme is also used 
to determine the local density, pressure, and entropy 
associated with each particle. With this setup, we 
have the benefit of being able to track thermodynamic 
quantities like in an SPH scheme but avoid the under 
mixing problem since the fluid equations are solved 
using the Eulerian scheme. 

Future work (Trac, Sills, & Pen 2003) will have 
higher resolution simulations. Collisions will be simu- 
lated in a box with 1024 x 1024 x 512 cells. Each parent 
star will have a radius of 192 grid cells and be assigned 
256 3 test particles. We will also be doing a detailed 
comparison between Eulerian and SPH simulations of 
stellar mergers. 

The self-gravitating hydro code used for the simu- 
lations is very memory friendly. For the 1024 2 x 512 
grid, 10 GB is required to store the hydro variables, 2 
GB for the potential, and less than 1 GB for the test 
particles. For every double timestep, approximately 
1000 floating point operations per grid cell is needed to 
carry out the TVD hydro calculations. The potential 
is computed once for every double step and this re- 



quires two FFTs. Since Eulerian codes are very mem- 
ory friendly, have low floating point counts, are eas- 
ily parallelized, and scale very well on shared-memory, 
multiple-processor machines, they can be used to run 
very high resolution simulations. 

8. Summary 

We have presented several numerical schemes for 
solving the linear advection equation and given the 
CFL stability conditions for each scheme. We have 
implemented the relaxing TVD scheme to solve the 
Euler system of conservation laws. The second-order 
accurate TVD scheme provides high resolution captur- 
ing of shocks, as can be seen in the Riemann shock test 
and the Sedov- Taylor blast wave test. The 1-D relax- 
ing TVD scheme can be easily generalized to higher 
dimensions using the dimensional splitting technique. 
A dimensionally split scheme can use longer time steps 
and is straightforward to implement in parallel. We 
have presented a sample astrophysical application. A 
3-D self-gravitating Eulerian hydro code is used to sim- 
ulate the formation of blue straggler stars through stel- 
lar mergers. We hope to have convinced the reader 
that Eulerian computational fluid dynamics is a pow- 
erful approach to simulating complex fluid flows be- 
cause it is simple, fast, and accurate. 

We thank Joachim Stadel and Norm Murray for 
comments and suggestions on the writing and editing 
of this paper. We also thank Alison Sills, Phil Arras, 
and Chris Matzner for discussions on stellar mergers. 
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A. 3-D Relaxing TVD code 



We provide a sample 3-D relaxing TVD code written in Fortran 90. The code is implemented using OpenMP 
directives to run in parallel on shared memory machines. The code is fast and memory friendly. The array u(a,i,j,k) 
stores the five conserved hydro quantities a = (p, pv x , pv y , pv z , e) for each cell (i,j,k) in the Cartesian cubical lattice 
with side length nc. For each sweep, we first call the subroutine timestep to determine the appropriate time step dt 
which satisfies the CFL condition. The updating of u by including the flux in the x direction is performed by the 
sweepx subroutine. The data array u is divided into 1-D array sections uld(a.i) which are operated on by the relaxing 
subroutine. The independent columns are distributed amongst multiple processors on a shared memory machine by 
the OpenMP directives. 

The relaxing TVD subroutine in this sample code is written for ease of readability and therefore, is not fully 
optimized. At the beginning of each partial step in the Runge-Kutta time integration scheme, the cell-averaged 
variables defined at grid cell centres are calculated by the averageflux subroutine. The fluxes at cell boundaries for 
the right-moving and left-moving waves are stored in fr and fl, respectively. We have implemented the minmod, 
superbee, and Van Leer flux limiters and the user of the code can easily switch between them. 

We have provided some initial conditions for the Sedov- Taylor blast wave test. The reader is encouraged to test the 
code and compare how the various flux limiters do at resolving strong shocks. This sample code does not implement 
the modified relaxing TVD scheme described at the end of §5.2, which has been found work very well with the Van 
Leer flux limiter but unstable with superbee for the 3-D Sedov Taylor test. We have found that the superbee limitcr 
is often unstable for 3-D fluid simulations. Please contact the authors regarding any questions on the implementation 
of the relaxing TVD algorithm. 

program main 
implicit none 

integer, parameter :: nc=64,hc=nc/2 
real, parameter :: gamma=5./3,cfl=0.9 

real, dimension(5,nc,nc,nc) :: u 

integer nsw,stopsim 
real t,tf,dt,E0,rmax 

t=0 
dt=0 
nsw=0 
stopsim=0 

E0=le5 
rmax=3*hc/4 

tf=sqrt((rmax/1.15)**5/E0) 

call sedovtaylor 
do 

call timestep 
call sweepx 
call sweepy 
call sweepz 
call sweepz 
call sweepy 
call sweepx 



18 



if (stopsim .eq. 1) exit 



call 


timestep 


call 


sweepz 


call 


sweepx 


call 


sweepy 


call 


sweepy 


call 


sweepx 


call 


sweepz 


if (stopsim .eq 


call 


timestep 


call 


sweepy 


call 


sweepz 


call 


sweepx 


call 


sweepx 


call 


sweepz 


call 


sweepy 


if (stopsim .eq 


enddo 





. 1) exit 



. 1) exit 



call outputresults 



contains 



subroutine sedovtaylor 
implicit none 
integer i,j,k 

do k=l,nc 
do j=l,nc 
do i=l,nc 
u(l,ij,k)=l 
u(2:4,i,j,k)=0 
u(5,i,j,k)=le-3 
enddo 
enddo 
enddo 

u(5,hc,hc,hc)=E0 
return 

end subroutine sedovtaylor 

subroutine outputresults 
implicit none 
integer i,j,k 
real r,x,y,z 

open(l,file='sedovtaylor.dat',recl=200) 
do k=l,nc 
z=k-hc 



do j=l,nc 
y=j-hc 
do i=l,nc 
x=i-hc 

r=sqrt(x**2+y**2+z**2) 
write(l,*) r,u(:,i,j,k) 
enddo 
enddo 
enddo 
close(l) 
return 

end subroutine outputresults 

subroutine timestep 
implicit none 
integer i,j,k 
real P.cs.cmax 
real v(3) 

cmax=le-5 

!$omp parallel do default(shared) private(i j,k,v,cs,P) reduction(max:cmax) 
do k=l,nc 
do j=l,nc 
do i=l,nc 

v=abs(u(2:4,i,j l k)/u(l,i,j,k)) 

P=max((gamma-l)*(u(5,i,j,k)-u(l,i,j,k)*sum(v**2)/2),0.) 
cs=sqrt(gamma*P/u(l,i,j,k)) 
cmax=max(cmax,maxval(v)+cs) 
enddo 
enddo 
enddo 

!$omp end parallel do 

dt=cfl/cmax 

if (t+2*dt .gt. tf) then 

dt=(tf-t)/2 

stopsim=l 
endif 

t=t+2*dt 
nsw=nsw+l 

write(*,"(a7,i3,a8,f7.5,a6,f8.5)") 'nsw = ',nsw,' dt = ',dt,' t = ',t 
return 

end subroutine timestep 

subroutine sweepx 
implicit none 
integer j,k 
real uld(5,nc) 
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!$omp parallel do default(shared) private(j,k,uld) 
do k=l,nc 
do j=l,nc 
uld=u(:,:,j,k) 
call relaxing(uld) 
u(:,:,j,k)=uld 
enddo 
enddo 

!$omp end parallel do 
return 
end subroutine sweepx 

subroutine sweepy 
implicit none 
integer i,k 
real uld(5,nc) 

!$omp parallel do default(shared) private(i,k,uld) 
do k=l,nc 
do i=l,nc 

uld((/l,3,2,4,5/),:)=u(:,i,:,k) 
call relaxing(uld) 
u(:,i,:,k)=uld((/l,3,2,4,5/),:) 
enddo 
enddo 

!$omp end parallel do 
return 
end subroutine sweepy 

subroutine sweepz 
implicit none 
integer i,j 
real uld(5,nc) 

!$omp parallel do default(shared) private(i,j,uld) 
do j=l,nc 
do i=l,nc 

uld((/l,4,3,2,5/),:)=u(:,i,j,:) 
call relaxing(uld) 
u(:,ij,:)=uld((/l,4,3,2,5/),:) 
enddo 
enddo 

!$omp end parallel do 
return 
end subroutine sweepz 
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subroutine relaxing(u) 
implicit none 
real, dimension(nc) :: c 
real, dimension(5,nc) :: u,ul,w,fu,fr,fl,dfl,dfr 

!! Do half step using first-order upwind scheme 
call averageflux(u,w,c) 
fr=(u*spread(c,l,5)+w)/2 
fl=cshift(u*spread(c,l,5)-w,l,2)/2 
fu=(fr-fl) 

ul=u-(fu-cshift(fu,-l,2))*dt/2 

!! Do full step using second-order TVD scheme 
call averageflux(ul,w,c) 

!! Right-moving waves 

fr=(ul*spread(c,l,5)+w)/2 

dfl=(fr-cshift(fr,-l,2))/2 

dfr=cshift(dfl,l,2) 

call vanleer(fr,dfl,dfr) 

kail minmod(fr,dfl,dfr) 

kail superbee(fr,dfl,dfr) 

!! Left-moving waves 

fl=cshift(ul*spread(c,l,5)-w,l,2)/2 

dfl=(cshift(fl,-l,2)-fl)/2 

dfr=cshift(dfl,l,2) 

call vanleer(fl,dfl,dfr) 

kail minmod(fl,dfl,dfr) 

kail superbee(fl,dfl,dfr) 

fu=(fr-fl) 

u=u-(fu-cshift(fu,-l,2))*dt 
return 

end subroutine relaxing 

subroutine averageflux(u,w,c) 
implicit none 
integer i 
real P,v 

real u(5,nc),w(5,nc),c(nc) 

!! Calculate cell-centered fluxes and freezing speed 
do i=l,nc 

v=u(2,i)/u(l,i) 

P=max((gamma-l)*(u(5,i)-sum(u(2:4,i)**2)/u(l,i)/2),0.) 

c(i)=abs(v)+max(sqrt(gamma*P/u(l,i)),le-5) 

w(l,i)=u(l,i)*v 

w(2,i)=(u(2,i)*v+P) 
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w(3,i)=u(3,i)*v 

w(4,i)=u(4,i)*v 

w(5,i)=(u(5,i)+P)*v 
enddo 
return 

end subroutine averageflux 

subroutine vanleer(f,a,b) 
implicit none 

real, dimension(5,nc) :: f,a,b,c 
c=a*b 

where (c .gt. 0) 

f=f+2*c/(a+b) 
endwhere 
return 
end subroutine vanleer 

subroutine minmod(f,a,b) 
implicit none 

real, dimension(nc) :: f,a,b 

f=f+(sign(l.,a)+sign(l.,b))*min(abs(a),abs(b))/2. 
return 

end subroutine minmod 



subroutine superbee(f,a,b) 
implicit none 

real, dimension(5,nc) :: f,a,b 

where (abs(a) .gt. abs(b)) 

f=f+(sign(l.,a)+sign(l.,b))*min(abs(a),abs(2*b))/2 
elsewhere 

f=f+(sign(l.,a)+sign(l.,b))*min(abs(2*a),abs(b))/2 
endwhere 
return 

end subroutine superbee 



end program main 



